Geometría Computacional

Gabriel De La Parra

Introducción

La propuesta inicial, es utilizar algoritmos de geometría en 2D para realizar procesamiento, análisis y extracción de features sobre distintas imágenes eléctricas.

Específicamente, se desea deseaba procesar distintos símbolos de diagramas eléctricos y obtener distintas métricas para cada uno. Dichas métricas pueden podrían haber servir para realizar el reconocimiento y clasificación de los mismos en otros escenarios.

Se buscó un enfoque distinto pero dentro del mismo campo, el manejo de imágenes con algoritmos geométricos para análisis y extracción de features.

En este trabajo se utilizarán los diagramas de Voronoi <informal> para lo que todo el mundo conoce que hacen, pero nadie hace.</informal> Se trabajará con la API de Google Maps para encontrar lugares de interés y sus diagramas de Voronoi con respecto a una dirección dada.

En una primera instancia se decidió ocupar C++ y CGAL para lo anterior, sin embargo se optó finalmente por trabajar en Python y SymPy.

CGAL es una librería muy poderosa, con diversas prestaciones y algoritmos en 2D y 3D. El autor considera que su instalación y uso es de mayor complejidad que otras opciones revisadas. De manera similar, al no tratarse de un escenario donde se requiere una alta capacidad computacional, no se justifica el uso de una librería tan compleja.

Dado lo anterior, se utilizará Python como lenguaje, por su simplicidad y su facilidad de integración con otras plataformas, como Jupyter Notebook, ipyWidgets y mpld3. Para los algoritmos geométricos se utilizará Sympy.

Sympy es una librería para matemática simbólica, sin embargo tiene varios módulos de geometría. Entre estos se encuentran implementaciones para:

  • Manejo de puntos, segmentos, polígonos
  • Cálculo de áreas e interesecciones
  • Cálculo de Cerradura Convexa (Convex Hull)
  • Triangulaciones
  • Diagramas de Voronoi
  • Cálculo del camino más corto

Desarrollo

El desarrollo de este trabajo se mostrará de forma incremental.

Se eligió Jupyter porque permite esta funcionalidad. De la misma manera, hay secciones interactivas, que permiten ajustar el resultado in-situ.

La siguiente celda demuestra el uso del código.

In [50]:
findThemAll("Beaucheff 851, Santiago", "sushi")

Imports

A continuación se presenta el listado de packages utilizados para esta experiencia.

In [1]:
"""Location"""
import ipywidgets as widgets
from IPython.display import clear_output
import googlemaps
gmaps = googlemaps.Client(key='AIzaSyDRGYgRRRgWbVEGzgmLXK4hf8jsxkBdpeI')

"""Voronoi"""
from io import BytesIO
from PIL import Image
from urllib import request
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (12,12)

"""Overlay"""
import warnings
warnings.filterwarnings('ignore')
from scipy.spatial import Voronoi, voronoi_plot_2d
import math
import numpy as np
from matplotlib.collections import LineCollection

"""Annotate"""
%matplotlib inline
from mpld3 import plugins
import mpld3

Geo Location

En primera instancia se definirá una clase GeoLocation para almacenar el nombre del lugar que se va a buscar y las coordenadas del mismo.

In [2]:
class GeoLocation:
    def __init__(self, name, lat, lng):
        self.name = name
        self.lat = lat
        self.lng = lng
    
    def __str__(self):
        return "{} ({},{})".format(self.name, self.lat, self.lng)
    
    def __repr__(self):
        return self.__str__()

Posteriormente se define un método para calcular el nombre y coordenadas según GoogleMaps para una búsqueda.

In [3]:
def getLocation(address):
    loc = gmaps.geocode(address)
    name = loc[0]["formatted_address"]
    lat = loc[0]["geometry"]["location"]["lat"]
    lng = loc[0]["geometry"]["location"]["lng"]
    return GeoLocation(name, lat, lng)
In [4]:
getLocation("Beaucheff 851, Santiago")
Out[4]:
Av. Beaucheff 851, Santiago, Región Metropolitana, Chile (-33.457624,-70.6642114)

Se presenta una búsqueda dinámica para pruebas.

In [5]:
def interactiveGetLocation(b):
    clear_output()
    address = locationTextBox.value
    location = getLocation(address)
    print(location)

locationButton = widgets.Button(description='Find')
locationButton.on_click(interactiveGetLocation)
locationTextBox = widgets.Textarea(value='Beaucheff 851, Santiago')
locationBox=widgets.VBox([locationTextBox, locationButton])

display(locationBox)

Visualizando el mapa

Ahora será turno de poder visualizar la dirección anterior. En este caso también se utilizará la API de Google Maps para obtener un mapa estático del lugar y colocar un marcador como referencia.

In [6]:
def getLocationImage(location, marker = True):
    url = "http://maps.googleapis.com/maps/api/staticmap?center={0},{1}&size=640x640&zoom=15".format(location.lat, location.lng)
    if marker:
        url += "&markers=color:red|{0},{1}".format(location.lat, location.lng)
    buffer = BytesIO(request.urlopen(url).read())
    return Image.open(buffer)

Adicionalmente, se separan los métodos de obtener la imagen de la visualización del mismo. Esto ya que más adelante ocuparemos esta imagen por si sola.

In [7]:
def printLocationMap(image):
    fig, ax = plt.subplots(figsize=(10,10))
    ax.imshow(image)
    plt.show()

Un ejemplo con lo anterior.

In [8]:
printLocationMap(getLocationImage(getLocation(locationTextBox.value), marker=True))

Buscando lugares

Una vez obtenida la ubicación, procederemos a encontrar los puntos de interés cercanos. Nuevamente definiremos un objeto para almacenar el nombre del lugar, su dirección y las coordenadas del mismo.

In [9]:
class NearbyPlace:
    def __init__(self, name, address, lat, lng):
        self.name = name
        self.address = address
        self.lat = lat
        self.lng = lng
    
    def __str__(self):
        return "{}\n{} ({},{})\n".format(self.name, self.address, self.lat, self.lng)
    
    def __repr__(self):
        return self.__str__()

Ahora un método para encontrar los lugares. Se utiliza nuevamente la API de Google Maps.

In [10]:
def searchPlaces(location, nearbyPlace):
    places = []
    local = gmaps.places_nearby(location=(location.lat,location.lng), radius=1000, keyword=nearbyPlace)
    for place in local["results"]:
        lat = place["geometry"]["location"]["lat"]
        lng = place["geometry"]["location"]["lng"]
        address = place["vicinity"]
        name = place["name"]
        places.append(NearbyPlace(name, address, lat, lng))
    return places

Un método interactivo para hacer pruebas. En este caso, con 'Restaurante'

In [11]:
def interactiveSearchPlaces(b):
    clear_output()
    nearbyPlace = searchPlaceTextBox.value
    places = searchPlaces(getLocation(locationTextBox.value), nearbyPlace)
    for place in places:
        print(place)

searchPlacesButton = widgets.Button(description='Find')
searchPlacesButton.on_click(interactiveSearchPlaces)
searchPlaceTextBox = widgets.Textarea(value='Restaurante')
searchBox=widgets.VBox([searchPlaceTextBox, searchPlacesButton])
display(searchBox)

Lugares en el mapa

Una vez obtenido los lugares de interés, se procederá a marcarlos y visualizarlos en el mapa.

In [12]:
def getPlacesImage(location, places, marker=True):    
    url = "http://maps.googleapis.com/maps/api/staticmap?center={0},{1}&size=640x640&zoom=15".format(location.lat, location.lng)
    if marker:
        url += "&markers=color:red|{0},{1}".format(location.lat, location.lng)
        for place in places:
            url += "&markers=color:blue|{0},{1}".format(place.lat,place.lng)
    buffer = BytesIO(request.urlopen(url).read())
    return Image.open(buffer)
In [13]:
location = getLocation(locationTextBox.value)
places = searchPlaces(location, searchPlaceTextBox.value)
printLocationMap(getPlacesImage(location, places))

Convirtiendo coordenadas

Vamos a pasar las coordenadas obtenidas a un plot. Para revisar como van quedando.

In [14]:
def getLocationPoints(location, places):
    pointsX = []
    pointsY = []
    
    for place in places:
        pointsX.append(place.lat-location.lat)
        pointsY.append(place.lng-location.lng)

    return pointsX, pointsY    
    
pointsX, pointsY = getLocationPoints(location, places)
plt.plot(pointsY, pointsX, 'ro', 0, 0, 'bs')
plt.show()

Voronoi de puntos

Una nota sobre los puntos anteriores. Si bien la distribución de los puntos es similar a la anterior, se puede apreciar que la escala está en latitud y longitud. Se deberá realizar una conversión posterior para sobreponerla sobre el mapa, sin embargo se verá más adelante.

Ahora se realizará y visualizará el diagrama de Voronoi sobre los puntos anteriores. Algo extraño con el algoritmo de la librería SymPy es que el diagrama toma los puntos en Y,X, por lo que se pasa como primer argumento la longitud y latitud.

In [15]:
def getPlacesPoints(pointsX, pointsY):
    points = []
    for i in range(len(pointsX)):
        points.append((pointsY[i], pointsX[i]))

    return points

pointsX, pointsY = getLocationPoints(location, places)
points = getPlacesPoints(pointsX, pointsY)
vor = Voronoi(points)
fig, ax = plt.subplots(figsize=(10, 10))
voronoi_plot_2d(vor, ax=ax)
plt.show()

Sobreponer Mapa y Voronoi

Como se mencionó anteriormente, los puntos de latitud y longitud están en una escala incorrecta.

Adicionalmente se deberá recentrarlos, ya que el cero en el mapa corresponde al centro y no a el borde de la imagen.

In [16]:
def relocatePlacesPoints(pointsX, pointsY):
    pointsX = [i*-28000+320 for i in pointsX]
    pointsY = [i*23000+320 for i in pointsY]
    return pointsX, pointsY
In [17]:
def drawVoronoiLines(points, ax):
    vor = Voronoi(points)
    line_colors = 'k'
    line_width = 1.0
    line_alpha = 1.0
    
    ax.set_autoscaley_on(False)
    ax.set_autoscalex_on(False)

    line_segments = []
    for simplex in vor.ridge_vertices:
        simplex = np.asarray(simplex)
        if np.all(simplex >= 0):
            line_segments.append([(x, y) for x, y in vor.vertices[simplex]])

    lc = LineCollection(line_segments,
                        colors=line_colors,
                        lw=line_width,
                        linestyle='solid')
    ax.add_collection(lc)
    ptp_bound = vor.points.ptp(axis=0)

    line_segments = []
    center = vor.points.mean(axis=0)
    for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices):
        simplex = np.asarray(simplex)
        if np.any(simplex < 0):
            i = simplex[simplex >= 0][0]  # finite end Voronoi vertex

            t = vor.points[pointidx[1]] - vor.points[pointidx[0]]  # tangent
            t /= np.linalg.norm(t)
            n = np.array([-t[1], t[0]])  # normal

            midpoint = vor.points[pointidx].mean(axis=0)
            direction = np.sign(np.dot(midpoint - center, n)) * n
            far_point = vor.vertices[i] + direction * ptp_bound.max()

            line_segments.append([(vor.vertices[i, 0], vor.vertices[i, 1]),
                                  (far_point[0], far_point[1])])

    lc = LineCollection(line_segments,
                        colors=line_colors,
                        lw=line_width,
                        linestyle='dashed')
    ax.add_collection(lc)
In [18]:
def overlayVoronoiMap(location, places):
    fig, ax = plt.subplots(figsize=(10,10))
    
    image = getLocationImage(location, marker=False)
    ax.imshow(image)

    pointsX, pointsY = getLocationPoints(location, places)
    pointsX, pointsY = relocatePlacesPoints(pointsX, pointsY)
    ax.plot(pointsY, pointsX, 'ro', 320, 320, 'bs')
    
    points = getPlacesPoints(pointsX, pointsY)
    drawVoronoiLines(points, ax)
       
    #Bonus:
    scatter = ax.scatter(x=pointsY, y=pointsX, alpha=0.1, s=1000)
    labels = ['<b><p>{}</p><p>{} ({},{})</p></b>'.format(i.name,i.address,i.lat,i.lng) for i in places]
    tooltip = mpld3.plugins.PointHTMLTooltip(scatter, labels=labels, css="b {color:black; font-size: 1em; display: block; text-shadow: 2px 0 0 #fff, -2px 0 0 #fff, 0 2px 0 #fff, 0 -2px 0 #fff, 1px 1px #fff, -1px -1px 0 #fff, 1px -1px 0 #fff, -1px 1px 0 #fff;}")
    mpld3.plugins.connect(fig, tooltip)
    mpld3.enable_notebook()
    
overlayVoronoiMap(location, places)
In [19]:
def findThemAll(whereYouAre, whatYouSearch):
    location = getLocation(whereYouAre)
    places = searchPlaces(location, whatYouSearch)
    overlayVoronoiMap(location, places)
In [51]:
findThemAll("Beaucheff 851, Santiago", "pizza")
In [37]:
findThemAll("Beaucheff 851, Santiago", "farmacia")
In [49]:
findThemAll("Beaucheff 851, Santiago", "cajero")